PREPRINT 



April 18, 2013 



Is there a hot spot of sea level rise acceleration along the 

mid- Atlantic United States? A Gaussian process decomposition 

of tide gauge records 

Robert E. Kopp 



o 

(N 

Oh 
< 

0^ 



Oh 
6 

O 



> 

o 

in 

^' 

o 
m 



X 



Recent research has suggested faster-than-global acceler- 
ation of sea level rise along the northeastern coast of the 
United States. To test this claim, I construct a Gaussian 
process model of tide gauge records, decomposing the tide 
gauge data into a smooth, long-term trend and red noise- 
type short-term variability (due to short-term ocean dynam- 
ics). This approach accommodates data gaps and allows de- 
composition of signals into local, regional, and global com- 
ponents. While tide gauge records do indicate a current 
faster-than-global increase in the rate of sea level rise in the 
mid-Atlantic region, this regional acceleration could reflect 
either the start of a long-term trend or short-term ocean 
dynamics. Current tide gauge data are insufficient to dis- 
criminate between these alternatives, and they will continue 
to be so for about two decades. 

1. Introduction 

Several recent papers report an increase in the rate of 
sea level rise (i.e., a sea level rise acceleration) along the 
mid- Atlantic coast of the United States that is greater than 
the global average [Sallenger et al., 2012; Boon, 2012; Ezer 
et al., 2013]. A plausible physical mechanism for producing 
such an acceleration exists: a slowing Gulf Stream should 
reduce the dynamic sea level gradient across the North At- 
lantic, causing a sea level decline in the central north At- 
lantic and a corresponding sea level rise in the northwestern 
north Atlantic. A secular sea level rise associated with this 
mechanism is one of the few points of agreement in global 
climate model projections of regional sea level change [Yin 
et al., 2009]. Ezer et al. [2013] suggested that changes in the 
Gulf Stream, reflected in the altimetry-observed alitudinal 
gradient across the Gulf Stream, are indeed correlated with 
a pronounced sea level rise in the eastern U.S. tide gauges 
since 2007. However, satellite altimetry of sea surface height 
dates back only a couple decades, so they were unable to in- 
vestigate the relationship between the Gulf Stream gradient 
and tide gauges before 1993. 

Regional and local sea level differ from mean global sea 
level (GSL) due to several factors. Glacial isostatic adjust- 
ment (GIA), which is driven by the viscous response of the 
mantle to the changes in ice mass loads since the Last Glacial 
Maximum, is causing the collapse of the proglacial forebulge 
in the mid- Atlantic region and the uplift of the Hudson Bay 
area. Land subsidence in the mid- Atlantic region leads to a 
sea level rise, partially mitigated by the geoid shift caused by 
the emplacement of mass underneath Hudson Bay [Peltier, 
2004] . GIA is a slow, long-term process; rates of GIA-driven 
sea level change are close to linear on the scale of the obser- 
vational record. East of the Fall Line, which passes close to 
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New York City, Baltimore, and Washington, bedrock is over- 
lain by the Mesozoic and Cenozoic sediments of the Coastal 
Plain, which can subside due to natural compaction and 
therefore experience a faster long-term rate of sea level rise. 
In the southern region of the Chesapeake Bay, high rates of 
subsidence may also be attributable to brecciation by the 
late Eocene Chesapeake Bay impact [Poag, 1997]. 

On shorter timescales, regional sea level anomalies can 
arise from land ice melt and ocean dynamics [e.g., Kopp 
et al., 2010]. Melting land ice leads to a slower or negative 
sea level rise near the meltwater source, and an enhanced 
sea level rise far from the source [Farrell and Clark, 1976; 
Mitrovica et al, 2001]. Greenland melt will thus produce 
less-than-global sea level rise on the Eastern seaboard of the 
United States, while Antarctic melt will produce greater- 
than-global sea level rise [Mitrovica et al., 2001, 2009]. Re- 
gional anomalies also arise from changes in ocean dynamic 
factors such as Gulf Stream strength, as previously noted; 
these factors can undergo both strong interannual variabil- 
ity and lower frequency variations. Local sea level anoma- 
lies can also arise from direct anthropogenic effects such as 
groundwater withdrawal and dredging. 

Here we develop a new spatio-temporal statistical method 
for analyzing tide gauge data, with the goal of partially 
disaggregating the sources of sea level rise. The method, 
based upon Gaussian process regression [Rasmussen and 
Williams, 2006] (alternatively known as spatio-temporal 
krieging) is well-suited for investigating possible regional ac- 
celerations for several reasons. First, it is fully probabilistic, 
so measurement and inferential uncertainties are propagated 
through the entire analysis. Second, it models sea level as 
a spatio-temporal field, naturally identifying regions with 
coherent sea level signals and appropriately sharing infor- 
mation among neighboring sites in the calculation of poste- 
rior sea level estimates. Third, being Bayesian in derivation 
(though empirical, rather than fully, Bayesian in implemen- 
tation), it copes naturally with the missing data that char- 
acterize tide gauge records. Fourth, it is non-parametric, 
and so does not force a functional form on the interpreta- 
tion of the tide gauge records (in contrast to the approaches 
of Sallenger et al. [2012] and Boon [2012]). It does, however, 
employ a parametric estimate of the prior covariance of sea 
level; this parameterized covariance allows easy separation 
of global, regional, and local signals, and of linear trends, 
smooth but non-linear variability, and red noise variability. 

Below, I demonstrate this method through application to 
tide gauge records from the eastern coast of North America 
and employ it to assess claims about a regional "hot spot" 
of sea level rise acceleration. 

2. Methods 

2.1. Statistical framework 

Sea level is a spatio-temporal field, which we model as a 
mean zero Gaussian process 



/(x,t)~eP(0,S) 



(1) 
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which we observe, through tide gauges, with independent 
but potentially heteroscedastic noise 



s(x,i) = /(x,i) + e(x,f). 



(2) 



The covariance S is the sum of nine terms, representing all 
combinations of three different temporal covariance func- 
tions and three different spatial covariance functions. For 
each term, the spatial and temporal covariances are separa- 
ble, though as the length scales of each spatial-temporal co- 
variance pair are different, they are not collectively so. The 
three temporal covariances have the forms of dot product 
(denoted by superscript '), rational quadratic ('), and expo- 
nential covariances ("). The dot product terms correspond 
to linear trends; the latter two terms are stationary, and rep- 
resent respectively smooth variations around the trend and 
high-frequency, red noise-type variability. The three spa- 
tial covariances represent globally uniform (g), regionally- 
specific (r), and site-specific components (i). 
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where C{r, v, 7) is a Matern covariance function of order 7, 
Axij is the angular distance between points i and j, Atij is 
the temporal distance ti — tj , and Sij is the Kronecker delta 
function. We interpret the regional, linear term as being due 
to GIA. I{i € CP,ji € CP) is an indicator function that is 
equal to 1 if sites i and j are both located on Coastal Plain 
sediments and if either site is located on bedrock; that 
is to say, for sites on bedrock we interpret all local linear 
trends as refiecting the regional GIA signal. 

We can condition upon observations s and hyperparam- 
eters 9 to estimate the posterior distribution of / at points 
(xt,t*), which we denote as f*: 



f.|s,e ~ AA(S..s(Ss + r,)-'s, 



(8) 



where T^ is the diagonal matrix of observation errors and 
Ss, S,, and S*,s represent respectively the covariance ma- 
trix for the points at which the observations s are taken, the 
covariance matrix for the desired points, and the covariance 
matrix between the desired points and the observed points. 
To isolate certain components of sea level in estimating f* 
(for example, to estimate only the component correspond- 
ing to the regional rational quadratic covariance function), 
we include in S, and St_s only the desired components. 
(See Fig. 1 for an example decomposition of the New York 
City (Battery) tide gauge record.) Note that this approach 
gives us a full posterior probability distribution, including 
the off-axis covariance terms. Accordingly, we can, through 
sampling, assess extreme value questions; for example, we 
can compare how a current rate of change compares to all 
rates of change during the last century. 



For the remainder of the paper, we use the term "sea 
level anomaly" to refer to the regional and local compo- 
nents of sea level {ar and (t„ terms), "regional sea level 
anomaly" to refer to the regional components of sea level ( 
(Jr terms), "smooth sea level anomaly" to refer to those com- 
ponents refiected in the temporally linear (S') and rational 
quadratic (S') terms, and "non-linear sea level anomaly" 
to refer to these components reflected in the temporally ra- 
tional quadratic and red noise-type terms (S"). A regional 
acceleration in sea level rise that is faster than the global 
average would accordingly be reflected by an accelerating 
rise in the smooth, regional sea level anomaly. 

2.2. Hyperparameters 

We denote the vector of spatial hyperparameters as 6x ~ 
{cr^,crJ.,crJj,7',cr«,cr?,cr^,7«,o-^,(7?,o-;j,7"} and the vector 
of temporal hyperparameters as 9t = {a,T'',T"'}. Rather 
than specifying priors over all of the hyperparameters, which 
would be the fully Bayesian approach, for computational ef- 
ficiency we adopt an empirical Bayesian approach, estimat- 
ing the values of the hyperparameters that maximizes the 
likelihood of the observations. For computational efficiency 
and interpretive ease, we adopt a staged approach. First, 
we optimize 9t individually for each tide gauge and then fix 
9t at the value of the median across tide gauges. To ensure 
that the selected hyperparameters separate red noise-type 
variability and smoother variations as intended, we bound 
r' < 1000 y and a > 0.05. Next, we fix the globally con- 
stant amplitudes {ag,a'^,ag} at the values that maximize 
the likelihood of the Church and White [2011] estimate of 
GSL. Third, we find the values of crj. and 7' that maximizes 
the likelihood of a Gaussian process, Matern covariance fit 
to the GIA rates of Peltier [2004] . Finally, we find the max- 
imum likelihood values of the remaining hyperparameters, 
using the median of the maximum likelihood parameters cal- 
culated with 20 data subsets. Each data subset consists 
of four mid-Atlantic tide gauges (Philadelphia, New York, 
Sandy Hook and Atlantic City) and twenty-five randomly 
selected 20-year blocks from other tide gauges. We employ 
high subset coverage in the mid- Atlantic because of the high 
spatial density of tide gauges in this region. Optimized hy- 
perparameters are shown in Table 1. 

2.3. Data 

We analyze estimates of mean annual sea level 
from the forty-seven eastern North America tide gauges 
archived by the Permanent Service for Mean Sea Level 
(http://www.psmsl.org/) with a record length exceeding 
thirty years. The gauges stretch from Daytona Beach, 
Florida, to St. John's, Newfoundland (Fig. SI). We in- 
directly incorporate data from other sites through the use 
of the Church and White [2011] estimate of GSL. Tide gauge 
observation errors are assumed to be ±6 mm (2a). 

3. Decomposition of tide gauge signals 

3.1. Long-term trends 

The linear trends projected at each site differ significantly 
from those projected using the ICE-5G VM2-90 glacio- 
isostatic adjustment model of Peltier [2004] (Figs. 2, S2). 
The discrepancy is particularly severe in Virginia, where the 
GIA projections lay outside the 95% confidence interval of 
the trend estimate. However, nowhere except between Mas- 
sachusetts and Pennsylvania and in South Carolina, north 
Florida, and parts of eastern maritime Canada do the ICE- 
5G projections fall within the 67% confidence interval of our 
inferred regional linear trends. Where data exists, our re- 
sults are generally in agreement (within uncertainty) with 
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the long-term rates of sea level rise estimated from multiple 
millennia of geological sea level proxies by Engelhart et al. 
[2009]. The notable exception to this agreement is on the 
Eastern Shore of the Chesapeake Bay, in Virginia, where 
the high linear trend seen in the tide gauge record is not 
matched in the geological record. 

The discrepancy between ICE-5G estimates and the ob- 
served long-term trends suggests either considerable error 
in the VM2-90 solid Earth parameters or the ICE-5G ice 
sheet history, or a major role for additional factors such as 
sediment compaction in the regions of discrepancy. Sedi- 
ment compaction may well be quite significant; the Virginia 
sites lay largely on Coastal Plain sediments in which com- 
paction could well play an important role. In addition to 
the Virginia sites, high linear rates of sea level rise relative 
to the projected regional trend are observed on the Coastal 
Plain of New Jersey at Sandy Hook and Atlantic City. The 
geographic spread of the high rate of sea level rise on the 
coastal plain throughout the mid-Atlantic region, not just 
in the vicinity of the Chesapeake Impact structure, suggests 
a dominant role for sediment compaction not related to the 
impact. Nevertheless, we do note considerable local variabil- 
ity among the four sites in the southern Chesapeake vicinity 
of the impact structure, with long term rates ranging from 
1.7 ± 0.8 (2cr) mm/y at Kipotopeke to 2.6 ± 0.6 mm/y at 
Sewells Point (Fig. S2). 

3.2. Non-linear component of regional sea level 

Fig. 3 shows the non-linear sea level anomaly. In Fig. 3d, 
four regionally coherent features stand out: significant rates 
of sea level anomaly rise along the entire seaboard in the 
1930s and 1940s; significant rates of sea level anomaly fall 
in the mid-Atlantic region in the 1970s, followed by slightly 
delayed fall to the north; and current significant rates of sea 
level anomaly rise in the mid- Atlantic and significant rates 
of sea level anomaly fall in the southeastern U.S. 

The sea level anomaly rise in the 1930s and 1940s is con- 
temporaneous with high rates of GSL rise, as can also be 
seen by examining the New York City tide gauge record 
(Fig. 1). The behavior of the Greenland ice sheet during 
this period is currently the subject of considerable disagree- 
ment [Gregory et al., 2013], with some modelers suggesting 
that the warm Northern Hemisphere temperatures of the 
1930s drove strong Greenland melt and others suggesting 
that it drove enhanced accumulation. The observed pattern 
of greater-than-global sea level rise off of North America dur- 
ing this interval is consistent with that expected from West 
Antarctic Ice Sheet melt and the opposite of what would be 
expected from the static sea level fingerprint of Greenland 
melt [Mitrovica et al., 2001, 2009]. However, the effects on 
the Gulf Stream of Greenland melt might be expected to 
operate in the opposite direction [Kopp et al, 2010]; it is 
therefore not possible from regional data alone to infer the 
role of Greenland in the 1930s and 1940s. A global analysis 
along the lines proposed by Hay et al. [2013] might succeed 
in this regard. 

In the mid- Atlantic region, the sea level anomaly fall in 
the 1970s and subsequent sea level anomaly rise in the 1990s 
and 2000s give rise to the "hot spot" of sea level rise accel- 
eration noted by other authors. We address the robustness 
and uniqueness of this feature below. 

4. A sea level rise hot spot? 

A weakening or migrating Gulf Stream will give rise to 
a reduction in the negative sea surface height gradient that 
runs north along the eastern North American coast [Yin 
et al., 2009]. In tide gauge records that use as a reference 



datum local sea level in a common year, this will appear as 
an increasing gradient between the northeast and the south- 
east. To evaluate the mid- Atlantic "hot spot," we therefore 
consider the difference in the non-linear regional sea level 
anomaly between New York City and Charleston (Fig. 4). 
The analysis removes the signals associated with global sea 
level change, with long-term linear changes due to effects 
such as GIA, and with purely local effects. 

Between 1990 and 2012, the smooth, non-linear regional 
sea level gradient between New York City and Charleston 
increased by 16±25 mm (an average rate of 0.7±1.1 mm/y); 
it is therefore very likely (probability ~ 90%) that the gra- 
dient has increased over this time period. The mid-Atlantic 
hot spot as such therefore does appear to be fairly robust. 
Its robustness does not, however, necessarily imply that the 
recent increases marks the start of a secular change in the 
Gulf Stream; it could reflect variability within the system. 

The greatest increase in the gradient over any 22-year pe- 
riod starting no earlier than 1900 and ending no later than 
1990 was 25 mm (95% range of 3-82 mm) (0.8 mm/y, range 
of 0.1-3.7 mm/y). It is likely (probability ~ 72%) that the 
rate of the 1990-2012 rise was exceeded at some point during 
the rest of the twentieth century. 

The magnitude of the gradient, referenced to the expected 
value in 1900 as a common datum, is currently 6 ± 38 mm, 
about 11 ± 26 mm below the maximum value attained be- 
tween 1900 and 1990. The current magnitude will need to 
increase by about 16 mm before it can be identified as likely 
(probability >~ 67%) unprecedented within the twentieth 
century and by about 28 mm before it can be identified as 
very likely unprecedented. At the average rate of increase 
of the last 22 years, this would take about 20 and 40 years, 
respectively. However, it is likely (probability ~ 80%) that 
the current rate of increase, estimated at 0.3 ± 2.2 mm/y, is 
less than the average over 1990-2012. 

Yin et al. [2009] project that, under the AlB scenario, a 
weakening Atlantic Meridional Overturning Circulation will 
give establish a ~ 150 mm dynamic sea level gradient be- 
tween New York and Miami during the century between 
1980-2000 and 2080-2100. The gradient between New York 
and Charleston will be similar. For such a gradient to be at- 
tained, another ~ 130 mm increase must occur over the next 
~ 80 years, on top of the ~ 16 mm that has occurred since 
1990. Achieving the gradient increase will require an accel- 
eration in the rate of increase of ~ 0.04 mm/y^. Comparing 
the mean rates of change over 1968-1990 to that over 1990- 
2012, we find an average acceleration of 0.05 ±0.08 mm/y^, 
which would likely be sufficient if sustained. It is about 
as likely as not (probability ~ 43%) that this acceleration 
was unprecedented in the rest of the twentieth century. If 
this acceleration were to continue for about two decades, the 
average rate of change would attain a likely unprecedented 
value of about 1.7 mm/y. 

While the current analysis is consistent with previous 
work identifying a recent shift to faster-than-global sea level 
rise in the mid-Atlantic region, neither the magnitude of 
the phenomenon, nor its rate of change, nor the accelera- 
tion of its rate of change appear to be beyond the bounds 
of twentieth-century variability. It is therefore premature 
to evaluate Sallenger et al. [2012] 's hypothesis that the cur- 
rent regionally high rates of sea level rise along the U.S. 
east coast represent the start of a long-term reorganization 
of the Gulf Stream, and it will take about two decades of 
additional observations before the sea level effects of such a 
reorganization can be identified in tide gauge records. 
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Figure 1. Example sea level decomposition, (a) Sea 
level at the Battery (blue), modeled underlying sea level 
field (red), and modeled underlying field without AR(1) 
component (green) . Church and White estimate of GSL 
shown for comparison (magneta). (b) The regionally and 
locally varying components of the sea level field at the 
battery, (c) The locally varying component of the sea 
level field, (d-f) corresponding 21-year running mean 
derivatives of the green curves in (a-c). The magneta 
curve in (d) is the corresponding GSL rate curve. Dashed 
and dotted lines show 67% and 95% confidence intervals, 
respectively. 
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Figure 2. (a) Map of mean estimate of long-term lin- 
ear sea level anomaly rate (i.e., change in sea level with 
respect to GSL). Dotted grey curves denotes the bound- 
aries of the Coastal Plain provinces \Fenneman and John- 
son, 1946]. (b) Regional linear sea level anomaly rates 
along the U.S. coast (blue; dashed=67% confidence inter- 
val, dotted=95% confidence interval), compared to ICE- 
5G projections of GIA-associated rates of sea level rise 
(red) [Peltier, 2004] and to the geological estimates of late 
Holocene sea level rise of Engelhart et al. [2009] (green; 
lines=l(T errors). 
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Figure 3. (a,b) Non-linear component of sea level 
anomaly and (c,d) smooth non-linear component of sea 
level anomaly rates. (b,d) show the regional component 
only, with (b) showing the expression and three spe- 
cific localities. In (a), heavy regions indicate space-time 
points where the standard deviation of the non-linear sea 
level anomaly estimate is < 20 mm. In (b), the green 
curves show the smooth component. Dashed lines denote 
67% confidence interval. In (c,d), heavy regions indicate 
space-time points where it is likely (probability > 67%) 
that the sign of the sea level anomaly rate component is 
correctly identified. 
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Figure 4. The difTerence in the regional sea level 
anomaly at New York City and Charleston, South Car- 
olina. Blue curves include red-noise type variability; 
green curves include only the smooth component, (a) 
Amplitude of the anomaly gradient, (b) rate of change. 
Dashed (dotted) lines denote 67% (95%) confidence in- 
tervals. 
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Supplementary text 

Optimized hyperparameters 

Optimized hyperparameters are shown in Table 1. The timescale of the red noise- type variabihty t" is 1 year. 
The rational quadratic parameters t'^ and a are fixed at their bounds, to facilitate interpretation; unconstrained 
optimization would lead the rational quadratic terms to incorporate some of the short-term variability accounted 
for b the red noise term. The length scales of the regional linear (7 ), rational quadratic (7"^) and red noise 
(7") terms are respectively 6.4, 4.2, and 7.6 degrees. No red noise term is needed for the global sea level curve 
(cr" — 0), which can be fully accommodated within its estimation error by smooth variability around a linear 
trend. At a regional scale, variability around the linear trend is accounted for in roughly equal proportion by 
smooth variations (a^ — 29.7 mm) and red noise (ct" = 28.8 mm). Smooth local variations are comparable in 
magnitude to regional variations {cti = 29.7 mm), while local red noise is smaller (tr" = 11.9 mm). 

Local non-linear sea level 

Examination of the non-linear local components highlights a few localities where the deviation of the non-linear 
components of local sea level from regional sea level are exceptionally high, as defined by the variance of the non- 
linear components. In particular, the five sites with at least half a century of data constituting top ten percent of 
variances are, in ranked order: St. John, New Brunswick; Philadelphia, Pennsylvania, Eastport, Maine; Portland, 
Maine; and Wilmington, North Carolina. Two of the sites are on the Gulf of Maine, and a third on the Bay of 
Fundy, and they are likely subject to high-amplitude, short wavelength dynamic variability. The Philadelphia 
tide gauge (Fig. S3) is located on the Delaware River, and the Wilmington tide gauge (Fig. S4) is located on the 
Cape Fear River. Both rivers have been dredged during the twentieth-century [Zervas, 2003]; the combination of 
dredging and natural riverine dynamics may account for the observed signals. 
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Figure SI. Tide gauge observations (blue) and the es- 
timated underlying sea level field (red; dotted=67% con- 
fidence interval). Green curves exclude red noise- type 
variability. Tide gauges presented from north to south. 
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Figure S2. Estimated linear component of sea level 
anomaly rate (blue — regional + local signal; green = 
regional signal only) compared to the ICE5G VM2-90 
GIA estimates. 
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Figure S3. Sea level decomposition for the tide gauge at Philadelphia. 



1900 1950 2000 





WILMINGTON 





a ji 


-50 


1 ; UF f 


-100 


1 ^Jfe/fl 


1 -150 


'ImJi'^ 


-200 


?;JttBl| 


-250 
-300 


/'^ 



Regional + Local 



Local 




1900 1950 2000 



1900 1950 2000 



1900 1950 2000 



* 3 

E 
£ 

S. 2 

? 1 

w 

-1 



a 




1900 1950 2000 




Figure S4. Sea level decomposition for the tide gauge at Wilmington. 



